GeneticsWithLongevity.R

#4 Генетика ~ Продолжительность жизни

rm(list=ls(all=TRUE))

#wd = getwd()
#wd = paste(wd, '/mtdna-mammalian-evolution/Body/2Derived',sep='')
#setwd(wd)

library(ggplot2)
data = read.table("../../Body/2Derived/CytB_with_ecology.csv", sep='\t', header = TRUE) 

###### КОРЕЛЯЦИИ

#(i) FractionOfSyn ~ -GenerLength; 

cor.test(x = data$GenerationLength_d, y = data$FractionOfSyn , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$FractionOfSyn, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$FractionOfSyn
## S = 10107033, p-value = 1.771e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.3786509

GenerationLength_FractionOfSyn_fit  <- cor.test(x = data$GenerationLength_d, y = data$FractionOfSyn)

ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfSyn ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfSyn , fill = Order)) + 
  geom_boxplot()

#(ii) FractionOfNonsyn ~ +GenerLength; 

cor.test(x = data$GenerationLength_d, y = data$FractionOfNonsyn , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$FractionOfNonsyn, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$FractionOfNonsyn
## S = 4555175, p-value = 1.771e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3786509

GenerationLength_FractionOfNonsyn_fit  <- cor.test(x = data$GenerationLength_d, y = data$FractionOfNonsyn)

ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfNonsyn ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfNonsyn , fill = Order)) + 
  geom_boxplot()

#(iii) SummOfAllGrantham ~ +GenerLength; 
cor.test(x = data$GenerationLength_d, y = data$SummOfAllGrantham , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$SummOfAllGrantham, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$SummOfAllGrantham
## S = 6351234, p-value = 0.01195
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1336593

GenerationLength_SummOfAllGrantham_fit  <- cor.test(x = data$GenerationLength_d, y = data$SummOfAllGrantham)

ggplot(data, aes(x = log(GenerationLength_d), y = SummOfAllGrantham ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = SummOfAllGrantham , fill = Order)) + 
  geom_boxplot()

#(iv) MedianOfAllNonsyn ~ +GenerLength 
cor.test(x = data$GenerationLength_d, y = data$MedianOfAllNonsyn , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$MedianOfAllNonsyn, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$MedianOfAllNonsyn
## S = 6992814, p-value = 0.3874
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04614444

GenerationLength_MedianOfAllNonsyn_fit  <- cor.test(x = data$GenerationLength_d, y = data$MedianOfAllNonsyn)

ggplot(data, aes(x = log(GenerationLength_d), y = MedianOfAllNonsyn ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y =MedianOfAllNonsyn , fill = Order)) + 
  geom_boxplot()

#(v) SummOfAllGrantham ~+SummOfAllNonsyn (проверка на вшивость)
cor.test(x = data$SummOfAllGrantham, y = data$MedianOfAllNonsyn, method = "spearm")

## Warning in cor.test.default(x = data$SummOfAllGrantham, y =
## data$MedianOfAllNonsyn, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$SummOfAllGrantham and data$MedianOfAllNonsyn
## S = 770600, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8948862

GenerationLength_MedianOfAllNonsyn_fit  <- cor.test(x = data$SummOfAllGrantham, y = data$MedianOfAllNonsyn)

ggplot(data, aes(x = log(SummOfAllGrantham), y = MedianOfAllNonsyn ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(SummOfAllGrantham), y = MedianOfAllNonsyn , fill = Order)) + 
  geom_boxplot()

## Warning: Removed 25 rows containing non-finite values (stat_boxplot).

#(vii) AverageGrantham ~ +GenerLength; 
cor.test(x = data$GenerationLength_d, y = data$AverageGrantham , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$AverageGrantham, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$AverageGrantham
## S = 5006904, p-value = 1.106e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3170328

GenerationLength_AverageGrantham_fit  <- cor.test(x = data$GenerationLength_d, y = data$AverageGrantham)

ggplot(data, aes(x = log(GenerationLength_d), y = AverageGrantham ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = AverageGrantham, fill = Order)) + 
  geom_boxplot()

#(viii) KnKs ~ +GenerLength; 
cor.test(x = data$GenerationLength_d, y = data$KnKs , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y = data$KnKs, : Cannot
## compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$KnKs
## S = 4588296, p-value = 3.592e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.374133

GenerationLength_KnKs_fit  <- cor.test(x = data$GenerationLength_d, y = data$KnKs)

ggplot(data, aes(x = log(GenerationLength_d), y = KnKs ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = KnKs, fill = Order)) + 
  geom_boxplot()

#(ix) DnDs ~ +GenerLength;
cor.test(x = data$GenerationLength_d, y = data$DnDs , method = "spearm")

## Warning in cor.test.default(x = data$GenerationLength_d, y = data$DnDs, : Cannot
## compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data$GenerationLength_d and data$DnDs
## S = 4135683, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4358718

GenerationLength_DnDs_fit  <- cor.test(x = data$GenerationLength_d, y = data$DnDs)

ggplot(data, aes(x = log(GenerationLength_d), y = DnDs ,col = factor(Order) ))+
  geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = DnDs, fill = Order)) + 
  geom_boxplot()

##################### ПИКИ

library(ape)
library(geiger)
library(caper)

## Loading required package: MASS

## Loading required package: mvtnorm

library(stringr)
library(phytools)

## Loading required package: maps

library(picante)

## Loading required package: vegan

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.5-7

## 
## Attaching package: 'vegan'

## The following object is masked from 'package:phytools':
## 
##     scores

## Loading required package: nlme

tree = read.tree("../../Body/1Raw/mammalia_species.nwk") 
data$Species = sub(" ", "_", data$Species, ignore.case = FALSE, perl = FALSE, fixed = FALSE, useBytes = FALSE)

row.names(data) = data$Species

nrow(data) # 344

## [1] 353

tree_w = treedata(tree, data[, c('Species', 'AverageGrantham', 'KnKs','DnDs' ,'GenerationLength_d')], 
                  sort=T, warnings=T)$phy

## Warning in treedata(tree, data[, c("Species", "AverageGrantham", "KnKs", : The following tips were not found in 'data' and were dropped from 'phy':
##  Abeomelomys_sevia
##  Abrocoma_bennettii
##  Abrocoma_cinerea
##  Abrothrix_andinus
##  Abrothrix_hershkovitzi
##  Abrothrix_illuteus
##  Abrothrix_jelskii
##  Abrothrix_lanosus
##  Abrothrix_markhami
##  Acerodon_celebensis
##  Acerodon_jubatus
##  Acinonyx_jubatus
##  Acomys_chudeaui
##  Acomys_cilicicus
##  Acomys_cineraceus
##  Acomys_ignitus
##  Acomys_johannis
##  Acomys_kempi
##  Acomys_minous
##  Acomys_nesiotes
##  Acomys_ngurui
##  Acomys_percivali
##  Acomys_russatus
##  Acomys_spinosissimus
##  Acomys_subspinosus
##  Acomys_wilsoni
##  Aconaemys_fuscus
##  Aconaemys_porteri
##  Aconaemys_sagei
##  Acrobates_pygmaeus
##  Addax_nasomaculatus
##  Aegialomys_xanthaeolus
##  Aepeomys_lugens
##  Aepyceros_melampus
##  Aepyprymnus_rufescens
##  Aeretes_melanopterus
##  Aeromys_tephromelas
##  Aethalops_aequalis
##  Aethalops_alecto
##  Aethomys_chrysophilus
##  Aethomys_ineptus
##  Aethomys_kaiseri
##  Ailurops_ursinus
##  Ailurus_fulgens
##  Akodon_aerosus
##  Akodon_affinis
##  Akodon_albiventer
##  Akodon_azarae
##  Akodon_baliolus
##  Akodon_boliviensis
##  Akodon_budini
##  Akodon_cursor
##  Akodon_dayi
##  Akodon_dolores
##  Akodon_fumeus
##  Akodon_iniscatus
##  Akodon_juninensis
##  Akodon_kofordi
##  Akodon_latebricola
##  Akodon_lindberghi
##  Akodon_lutescens
##  Akodon_mimus
##  Akodon_molinae
##  Akodon_mollis
##  Akodon_mystax
##  Akodon_orophilus
##  Akodon_paranaensis
##  Akodon_reigi
##  Akodon_serrensis
##  Akodon_siberiae
##  Akodon_simulator
##  Akodon_subfuscus
##  Akodon_sylvanus
##  Akodon_toba
##  Akodon_torques
##  Akodon_varius
##  Alcelaphus_buselaphus
##  Alcelaphus_caama
##  Alcelaphus_lichtensteini
##  Alces_alces
##  Alces_americanus
##  Alionycteris_paucidentata
##  Allactaga_balikunica
##  Allactaga_bullata
##  Allactaga_elater
##  Allactaga_euphratica
##  Allactaga_firouzi
##  Allactaga_hotsoni
##  Allactaga_major
##  Allactaga_sibirica
##  Allactaga_sp._Kazakhstan
##  Allactaga_toussi
##  Allactaga_williamsi
##  Allactodipus_bobrinskii
##  Allenopithecus_nigroviridis
##  Allocebus_trichotis
##  Allocricetulus_curtatus
##  Allocricetulus_eversmanni
##  Alouatta_palliata
##  Alouatta_pigra
##  Alouatta_sara
##  Alouatta_seniculus
##  Alticola_albicaudus
##  Alticola_argentatus
##  Alticola_barakshin
##  Alticola_macrotis
##  Alticola_semicanus
##  Alticola_strelzowi
##  Alticola_tuvinicus
##  Amblysomus_corriae
##  Amblysomus_hottentotus
##  Amblysomus_marleyi
##  Amblysomus_septentrionalis
##  Ametrida_centurio
##  Ammospermophilus_harrisii
##  Ammospermophilus_insularis
##  Ammospermophilus_interpres
##  Ammospermophilus_leucurus
##  Ammotragus_lervia
##  Amphinectomys_savamis
##  Anathana_ellioti
##  Andalgalomys_olrogi
##  Andalgalomys_pearsoni
##  Andalgalomys_roigi
##  Andinomys_edax
##  Anisomys_imitator
##  Anomalurus_beecrofti
##  Anoura_caudifer
##  Anoura_cultrata
##  Anoura_geoffroyi
##  Anoura_latidens
##  Anourosorex_yamashinai
##  Antechinomys_laniger
##  Antechinus_agilis
##  Antechinus_bellus
##  Antechinus_flavipes
##  Antechinus_godmani
##  Antechinus_leo
##  Antechinus_minimus
##  Antechinus_stuartii
##  Antechinus_swainsonii
##  Anthops_ornatus
##  Antidorcas_marsupialis
##  Antilocapra_americana
##  Antilope_cervicapra
##  Antrozous_dubiaquercus
##  Antrozous_pallidus
##  Aonyx_capensis
##  Aonyx_cinerea
##  Aotus_azarai
##  Aotus_brumbacki
##  Aotus_griseimembra
##  Aotus_lemurinus
##  Aotus_nancymaae
##  Aotus_nigriceps
##  Aotus_trivirgatus
##  Aotus_vociferans
##  Aplodontia_rufa
##  Apodemus_alpicola
##  Apodemus_epimelas
##  Apodemus_gurkha
##  Apodemus_hermonensis
##  Apodemus_hyrcanicus
##  Apodemus_iconicus
##  Apodemus_ilex
##  Apodemus_mystacinus
##  Apodemus_pallipes
##  Apodemus_ponticus
##  Apodemus_semotus
##  Apodemus_uralensis
##  Apodemus_wardi
##  Apodemus_witherbyi
##  Apomys_abrae
##  Apomys_camiguinensis
##  Apomys_datae
##  Apomys_gracilirostris
##  Apomys_hylocoetes
##  Apomys_insignis
##  Apomys_iridensis
##  Apomys_lubangensis
##  Apomys_microdon
##  Apomys_musculus
##  Apomys_sacobianus
##  Apomys_sp._SJS-2010a
##  Apomys_sp._SJS-2010b
##  Apomys_sp._SJS-2010c
##  Apomys_sp._SJS-2010d
##  Apomys_sp._SJS-2010e
##  Apomys_sp._SJS-2010f
##  Apomys_sp._SJS-2010g
##  Aproteles_bulmerae
##  Arborimus_albipes
##  Arborimus_longicaudus
##  Arborimus_pomo
##  Archboldomys_luzonensis
##  Archboldomys_sp._SAJ-2012
##  Arctictis_binturong
##  Arctocebus_aureus
##  Arctocebus_calabarensis
##  Arctocephalus_australis
##  Arctocephalus_galapagoensis
##  Arctocephalus_gazella
##  Arctocephalus_philippii
##  Arctocephalus_pusillus
##  Arctocephalus_townsendi
##  Arctocephalus_tropicalis
##  Arctodus_simus
##  Arctogalidia_trivirgata
##  Arctonyx_collaris
##  Ardops_nichollsi
##  Arielulus_aureocollaris
##  Arielulus_circumdatus
##  Arielulus_cuprosus
##  Ariteus_flavescens
##  Artibeus_amplus
##  Artibeus_anderseni
##  Artibeus_aztecus
##  Artibeus_bogotensis
##  Artibeus_concolor
##  Artibeus_fraterculus
##  Artibeus_hartii
##  Artibeus_hirsutus
##  Artibeus_inopinatus
##  Artibeus_intermedius
##  Artibeus_schwartzi
##  Artibeus_toltecus
##  Arvicanthis_abyssinicus
##  Arvicanthis_ansorgei
##  Arvicanthis_nairobae
##  Arvicanthis_neumanni
##  Arvicanthis_niloticus
##  Arvicanthis_somalicus
##  Arvicola_amphibius
##  Arvicola_sapidus
##  Arvicola_scherman
##  Asellia_tridens
##  Aselliscus_stoliczkanus
##  Aselliscus_tricuspidatus
##  Atelerix_albiventris
##  Atelerix_algirus
##  Ateles_belzebuth
##  Ateles_fusciceps
##  Ateles_geoffroyi
##  Ateles_hybridus
##  Ateles_paniscus
##  Atelocynus_microtis
##  Atherurus_africanus
##  Atherurus_macrourus
##  Atilax_paludinosus
##  Atlantoxerus_getulus
##  Auliscomys_boliviensis
##  Auliscomys_micropus
##  Auliscomys_pictus
##  Auliscomys_sublimis
##  Avahi_cleesei
##  Avahi_occidentalis
##  Avahi_peyrierasi
##  Avahi_unicolor
##  Axis_axis
##  Axis_kuhlii
##  Axis_porcinus
##  Babyrousa_babyrussa
##  Baiomys_musculus
##  Baiomys_taylori
##  Balaenoptera_acutorostrata
##  Balaenoptera_bonaerensis
##  Balaenoptera_borealis
##  Balaenoptera_brydei
##  Balaenoptera_edeni
##  Balaenoptera_musculus
##  Balaenoptera_omurai
##  Balantiopteryx_infusca
##  Balantiopteryx_io
##  Balantiopteryx_plicata
##  Balionycteris_maculata
##  Bandicota_bengalensis
##  Bandicota_savilei
##  Barbastella_barbastellus
##  Barbastella_beijingensis
##  Barbastella_darjelingensis
##  Barbastella_leucomelas
##  Bassaricyon_alleni
##  Bassaricyon_beddardi
##  Bassaricyon_gabbii
##  Bassariscus_astutus
##  Bassariscus_sumichrasti
##  Bathyergus_janetta
##  Bathyergus_suillus
##  Batomys_granti
##  Batomys_salomonseni
##  Bdeogale_crassicauda
##  Bdeogale_nigripes
##  Beamys_hindei
##  Beatragus_hunteri
##  Belomys_pearsonii
##  Berardius_arnuxii
##  Berardius_bairdii
##  Berylmys_berdmorei
##  Berylmys_bowersi
##  Bettongia_gaimardi
##  Bettongia_lesueur
##  Bettongia_penicillata
##  Bettongia_tropica
##  Bibimys_chacoensis
##  Bibimys_labiosus
##  Blanfordimys_afghanus
##  Blanfordimys_bucharensis
##  Blarina_carolinensis
##  Blarina_hylophaga
##  Blarinomys_breviceps
##  Blastocerus_dichotomus
##  Bolomys_amoenus
##  Bolomys_lasiurus
##  Bolomys_temchuki
##  Bolomys_urichi
##  Boneia_bidens
##  Bos_frontalis
##  Bos_gaurus
##  Bos_grunniens
##  Bos_indicus
##  Bos_javanicus
##  Bos_mutus
##  Bos_primigenius
##  Bos_sauveli
##  Bos_taurus
##  Boselaphus_tragocamelus
##  Brachiones_przewalskii
##  Brachylagus_idahoensis
##  Brachyphylla_nana
##  Brachyphylla_pumila
##  Brachytarsomys_albicauda
##  Brachyteles_arachnoides
##  Brachyteles_hypoxanthus
##  Brachyuromys_betsileoensis
##  Brachyuromys_ramirohitra
##  Bradypus_pygmaeus
##  Bradypus_torquatus
##  Bradypus_tridactylus
##  Bradypus_variegatus
##  Brucepattersonius_igniventris
##  Brucepattersonius_soricinus
##  Bubalus_bubalis
##  Bubalus_carabanensis
##  Bubalus_depressicornis
##  Bubalus_mindorensis
##  Bubalus_quarlesi
##  Budorcas_taxicolor
##  Bullimus_bagobus
##  Bullimus_gamay
##  Bullimus_luzonicus
##  Bunolagus_monticularis
##  Bunomys_andrewsi
##  Bunomys_chrysocomus
##  Burramys_parvus
##  Cabassous_centralis
##  Cabassous_chacoensis
##  Cabassous_tatouay
##  Cabassous_unicinctus
##  Cacajao_ayresi
##  Cacajao_hosomi
##  Cacajao_melanocephalus
##  Caenolestes_fuliginosus
##  Calcochloris_obtusirostris
##  Callicebus_coimbrai
##  Callicebus_dubius
##  Callicebus_nigrifrons
##  Callicebus_personatus
##  Callimico_goeldii
##  Callithrix_argentata
##  Callithrix_aurita
##  Callithrix_emiliae
##  Callithrix_geoffroyi
##  Callithrix_humeralifera
##  Callithrix_humilis
##  Callithrix_jacchus
##  Callithrix_kuhlii
##  Callithrix_mauesi
##  Callithrix_melanura
##  Callithrix_penicillata
##  Callithrix_pygmaea
##  Callithrix_rondoni
##  Callithrix_saterei
##  Callorhinus_ursinus
##  Callosciurus_caniceps
##  Callosciurus_finlaysonii
##  Callosciurus_inornatus
##  Callosciurus_nigrovittatus
##  Callosciurus_notatus
##  Callosciurus_orestes
##  Callosciurus_prevostii
##  Callospermophilus_lateralis
##  Callospermophilus_madrensis
##  Callospermophilus_saturatus
##  Calomys_callidus
##  Calomys_fecundus
##  Calomys_hummelincki
##  Calomys_laucha
##  Calomys_lepidus
##  Calomys_musculinus
##  Calomys_sorellus
##  Calomys_sp.
##  Calomys_sp._CEG40
##  Calomys_tener
##  Calomys_tocantinsi
##  Calomys_venustus
##  Calomys

## Warning in treedata(tree, data[, c("Species", "AverageGrantham", "KnKs", : The following tips were not found in 'phy' and were dropped from 'data':
##  Crocidura_monticola
##  Eidolon_dupreanum
##  Epomops_buettikoferi
##  Sorex_cylindricauda

data_tree = as.data.frame(treedata(tree_w, data[, c('Species', 'AverageGrantham', 'KnKs','DnDs' ,'GenerationLength_d')], 
                                   sort=T, warnings=T)$data)

## Warning in treedata(tree_w, data[, c("Species", "AverageGrantham", "KnKs", : The following tips were not found in 'phy' and were dropped from 'data':
##  Crocidura_monticola
##  Eidolon_dupreanum
##  Epomops_buettikoferi
##  Sorex_cylindricauda

setdiff(data$Species, data_tree$Species)# 4 вида отвалились при присоединении дерева

## [1] "Crocidura_monticola"  "Eidolon_dupreanum"    "Epomops_buettikoferi"
## [4] "Sorex_cylindricauda"

data_tree$AverageGrantham = as.numeric(as.character(data_tree$AverageGrantham))
data_tree$GenerationLength_d = as.numeric(as.character(data_tree$GenerationLength_d))
data_tree$KnKs = as.numeric(as.character(data_tree$KnKs))

#GenLen & Grantham
cor.test(pic(data_tree$AverageGrantham, tree_w), pic(log2(data_tree$GenerationLength_d), tree_w), method = 'spearman')

## Warning in cor.test.default(pic(data_tree$AverageGrantham, tree_w),
## pic(log2(data_tree$GenerationLength_d), : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  pic(data_tree$AverageGrantham, tree_w) and pic(log2(data_tree$GenerationLength_d), tree_w)
## S = 6609013, p-value = 0.2717
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.05907785

cor.test(data_tree$AverageGrantham,data_tree$GenerationLength_d, method = 'spearman')

## Warning in cor.test.default(data_tree$AverageGrantham,
## data_tree$GenerationLength_d, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data_tree$AverageGrantham and data_tree$GenerationLength_d
## S = 4785484, p-value = 5.302e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3245325

#GenLen & DnDs
cor.test(pic(data_tree$DnDs, tree_w), pic(log2(data_tree$GenerationLength_d), tree_w), method = 'spearman')

## Warning in cor.test.default(pic(data_tree$DnDs, tree_w),
## pic(log2(data_tree$GenerationLength_d), : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  pic(data_tree$DnDs, tree_w) and pic(log2(data_tree$GenerationLength_d), tree_w)
## S = 6113720, p-value = 0.01556
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1295925

cor.test(as.numeric(data_tree$DnDs),data_tree$GenerationLength_d, method = 'spearman')

## Warning in cor.test.default(as.numeric(data_tree$DnDs),
## data_tree$GenerationLength_d, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  as.numeric(data_tree$DnDs) and data_tree$GenerationLength_d
## S = 3984711, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4375611

#GenLen & KnKs
cor.test(pic(data_tree$KnKs, tree_w), pic(log2(data_tree$GenerationLength_d), tree_w), method = 'spearman')

## Warning in cor.test.default(pic(data_tree$KnKs, tree_w),
## pic(log2(data_tree$GenerationLength_d), : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  pic(data_tree$KnKs, tree_w) and pic(log2(data_tree$GenerationLength_d), tree_w)
## S = 6580489, p-value = 0.2401
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.06313875

cor.test(data_tree$KnKs,data_tree$GenerationLength_d, method = 'spearman')

## Warning in cor.test.default(data_tree$KnKs, data_tree$GenerationLength_d, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  data_tree$KnKs and data_tree$GenerationLength_d
## S = 4441184, p-value = 5.695e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3731303

######### МАНН-УИТНИ
Data = data
quantile(Data$GenerationLength_d, seq(from=0, to=1, by=0.25),na.rm = TRUE)  

##        0%       25%       50%       75%      100% 
##   202.795   628.295  1585.196  2726.536 18980.000

boxplot(Data[Data$GenerationLength_d < 637.3329,]$AverageGrantham, +
          Data[Data$GenerationLength_d < 1825.0000 & Data$GenerationLength_d > 637.3329,]$AverageGrantham, +
          Data[Data$GenerationLength_d < 2869.6933 & Data$GenerationLength_d > 1825.0000,]$AverageGrantham, +
          Data[Data$GenerationLength_d > 2869.6933,]$AverageGrantham, names = c('Q1','Q2', 'Q3','Q4'), outline = FALSE, xlab = "GenLenght", ylab = "AverageGrantham",notch = TRUE)

boxplot(Data[Data$GenerationLength_d < 637.3329,]$KnKs, +
          Data[Data$GenerationLength_d < 1825.0000 & Data$GenerationLength_d > 637.3329,]$KnKs, +
          Data[Data$GenerationLength_d < 2869.6933 & Data$GenerationLength_d > 1825.0000,]$KnKs, +
          Data[Data$GenerationLength_d > 2869.6933,]$KnKs, names = c('Q1','Q2', 'Q3','Q4'), outline = FALSE, xlab = "GenLenght", ylab = "KnKs",notch = TRUE)

##### boxplots by quartiles
boxplot(data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
        data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham,
        data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$AverageGrantham,
        data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$AverageGrantham,
        names=c('1','2','3','4'), outline = FALSE, notch = TRUE)

boxplot(data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$DnDs,
        data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$DnDs,
        data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$DnDs,
        data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$DnDs,
        names=c('1','2','3','4'), outline = FALSE, notch = TRUE)

wilcox.test(data$AverageGrantham,data$GenerationLength_d)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data$AverageGrantham and data$GenerationLength_d
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(data$DnDs,data$GenerationLength_d)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data$DnDs and data$GenerationLength_d
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(data$KnKs,data$GenerationLength_d)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data$KnKs and data$GenerationLength_d
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(
  data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data[data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.25), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.25) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.5), ]$AverageGrantham
## W = 3011.5, p-value = 0.007993
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(
  data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$AverageGrantham)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data[data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.25), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.5) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 2848, p-value = 0.001206
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(
  data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$AverageGrantham)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data[data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.25), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 2153, p-value = 1.449e-07
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham,
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$AverageGrantham)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.25) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.5), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.5) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 3478.5, p-value = 0.2975
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham,
  data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$AverageGrantham)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.25) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.5), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 2621, p-value = 0.0003169
## alternative hypothesis: true location shift is not equal to 0

hist(data$GenerationLength_d,breaks=10)

######### ЛЯМБДА

GenLen <- data_tree$GenerationLength_d
names(GenLen) <- rownames(data_tree)

Grantham <- data_tree$AverageGrantham
names(Grantham) <- rownames(data_tree)

DnDs <- as.matrix(as.numeric(data_tree$DnDs))
rownames(DnDs) <- rownames(data_tree)

phylosig(tree_w, GenLen, method = "lambda", test = TRUE) 

## 
## Phylogenetic signal lambda : 0.993386 
## logL(lambda) : -2845.39 
## LR(lambda=0) : 630.524 
## P-value (based on LR test) : 3.84386e-139

phylosig(tree_w, Grantham, method = "lambda", test = TRUE)

## 
## Phylogenetic signal lambda : 0.196536 
## logL(lambda) : -1468.53 
## LR(lambda=0) : 17.9445 
## P-value (based on LR test) : 2.27441e-05

phylosig(tree_w, DnDs, method = "lambda", test = TRUE)

## 
## Phylogenetic signal lambda : 0.465285 
## logL(lambda) : 554.683 
## LR(lambda=0) : 79.2683 
## P-value (based on LR test) : 5.42225e-19

################### PGLS

MutComp = comparative.data(tree_w, data, Species, vcv=TRUE)

summary(pgls(scale(AverageGrantham) ~ log2(GenerationLength_d), MutComp, lambda="ML"))

## 
## Call:
## pgls(formula = scale(AverageGrantham) ~ log2(GenerationLength_d), 
##     data = MutComp, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.248868 -0.046084 -0.001791  0.044802  0.212354 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.000
##    lower bound : 0.000, p = 1    
##    upper bound : 1.000, p = < 2.22e-16
##    95.0% CI   : (NA, 0.293)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                           Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)              -2.089236   0.425581 -4.9091 1.411e-06 ***
## log2(GenerationLength_d)  0.199745   0.040386  4.9459 1.183e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07277 on 347 degrees of freedom
## Multiple R-squared: 0.06585, Adjusted R-squared: 0.06316 
## F-statistic: 24.46 on 1 and 347 DF,  p-value: 1.183e-06

summary(pgls(scale(KnKs) ~ log2(GenerationLength_d), MutComp, lambda="ML"))

## 
## Call:
## pgls(formula = scale(KnKs) ~ log2(GenerationLength_d), data = MutComp, 
##     lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32732 -0.05094 -0.00464  0.04000  0.31108 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.334
##    lower bound : 0.000, p = 0.00070586
##    upper bound : 1.000, p = < 2.22e-16
##    95.0% CI   : (0.086, 0.610)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                           Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)              -2.708434   0.870111 -3.1127 0.0020075 ** 
## log2(GenerationLength_d)  0.264797   0.075496  3.5074 0.0005119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07904 on 347 degrees of freedom
## Multiple R-squared: 0.03424, Adjusted R-squared: 0.03146 
## F-statistic:  12.3 on 1 and 347 DF,  p-value: 0.0005119

summary(pgls(scale(DnDs) ~ log2(GenerationLength_d), MutComp, lambda="ML"))

## 
## Call:
## pgls(formula = scale(DnDs) ~ log2(GenerationLength_d), data = MutComp, 
##     lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29521 -0.03523  0.00188  0.03274  0.31114 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.331
##    lower bound : 0.000, p = 3.17e-12
##    upper bound : 1.000, p = < 2.22e-16
##    95.0% CI   : (0.171, 0.534)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                           Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)              -3.109900   0.823284 -3.7774 0.0001864 ***
## log2(GenerationLength_d)  0.282557   0.071471  3.9535 9.337e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07493 on 347 degrees of freedom
## Multiple R-squared: 0.0431,  Adjusted R-squared: 0.04034 
## F-statistic: 15.63 on 1 and 347 DF,  p-value: 9.337e-05

######################## ДЕЛЕНИЕ ПО СЕМЕЙСТВАМ

OrderFreq = data.frame(table(data$Order));
FrequentOrders = OrderFreq[OrderFreq$Freq >= 3,]$Var1; length(FrequentOrders) # 3!!!

## [1] 8

OrderFreq = OrderFreq[OrderFreq$Var1 %in% FrequentOrders,]
names(OrderFreq)=c('Order','NumberOfSpecies')

Mammalia = data[data$Order %in% FrequentOrders,]
agg = aggregate(list(Mammalia$AverageGrantham,Mammalia$GenerationLength_d,Mammalia$DnDs,Mammalia$KnKs), by = list(Mammalia$Order), FUN = median)
names(agg) = c('Order','AverageGrantham','GenerationLength_d','DnDs','KnKs')

test1 = cor.test(agg$AverageGrantham,agg$GenerationLength_d,method = 'spearman') ### 0.008503, Rho = 0.8434347 !!!#не работает с пропусками в продолжительности жизни

## Warning in cor.test.default(agg$AverageGrantham, agg$GenerationLength_d, :
## Cannot compute exact p-value with ties

plot(agg$GenerationLength_d,agg$AverageGrantham)

test2 = cor.test(agg$DnDs,agg$GenerationLength_d,method = 'spearman') ### 0.01071, Rho = 0.8571429 !!!
plot(agg$DnDs,agg$AverageGrantham)

test3 = cor.test(agg$KnKs,agg$GenerationLength_d,method = 'spearman') ### 0.004563, Rho = 0.9047619 !!!
plot(agg$KnKs,agg$AverageGrantham)

OrderGrantham_GenLen = c(test1$p.value,test1$estimate)
OrderDnDs_GenLen = c(test2$p.value,test2$estimate)
OrderKnKs_GenLen = c(test3$p.value,test3$estimate)

OrderFinal = data.frame(OrderGrantham_GenLen,OrderDnDs_GenLen,OrderKnKs_GenLen); rownames(OrderFinal) = c('p-value','rho')

FamFreq = data.frame(table(data$family));
FrequentFamilies = FamFreq[FamFreq$Freq >= 3,]$Var1; length(FrequentFamilies) # 3!!!

## [1] 32

FamFreq = FamFreq[FamFreq$Var1 %in% FrequentFamilies,]
names(FamFreq)=c('Family','NumberOfSpecies')

# Family
Mammalia1 = data[data$family %in% FrequentFamilies,]
agg1 = aggregate(list(Mammalia1$AverageGrantham,Mammalia1$GenerationLength_d,Mammalia1$DnDs,Mammalia1$KnKs), by = list(Mammalia1$family), FUN = median)
names(agg1) = c('Family','AverageGrantham','GenerationLength_d','DnDs','KnKs')

test4 = cor.test(agg1$AverageGrantham,agg1$GenerationLength_d,method = 'spearman') ### 0.00281, Rho = 0.5108606 !!!
plot(agg1$GenerationLength_d,agg1$AverageGrantham)

test5 = cor.test(agg1$DnDs,agg1$GenerationLength_d,method = 'spearman') ### 0.003045, Rho = 0.5128299 !!!
plot(agg1$DnDs,agg1$AverageGrantham)

test6 = cor.test(agg1$KnKs,agg1$GenerationLength_d,method = 'spearman') ### 0.01945, Rho = 0.4109991  !!!

## Warning in cor.test.default(agg1$KnKs, agg1$GenerationLength_d, method =
## "spearman"): Cannot compute exact p-value with ties

plot(agg1$KnKs,agg1$AverageGrantham)

FamilyGrantham_GenLen = c(test4$p.value,test4$estimate)
FamilyDnDs_GenLen = c(test5$p.value,test5$estimate)
FamilyKnKs_GenLen = c(test6$p.value,test6$estimate)

FamilyFinal = data.frame(FamilyGrantham_GenLen,FamilyDnDs_GenLen,FamilyKnKs_GenLen); rownames(FamilyFinal) = c('p-value','rho')

agg = merge(agg,OrderFreq)
agg = agg[order(agg$GenerationLength_d),]##"Eulipotyphla","Rodentia","Didelphimorphia","Lagomorpha","Chiroptera","Carnivora","Cetartiodactyla","Primates"
agg$AverageGrantham = round(agg$AverageGrantham,2)
agg$GenerationLength_d = round(agg$GenerationLength_d,0)

#################### Боксплоты по отрядам и семействам

Order=c("Eulipotyphla","Rodentia","Didelphimorphia","Lagomorpha","Chiroptera","Carnivora","Cetartiodactyla","Cetartiodactyla")#    
median_DnDs = c(median(Data[Data$Order == 'Eulipotyphla',]$DnDs),median(Data[Data$Order == 'Rodentia',]$DnDs),median(Data[Data$Order == 'Didelphimorphia',]$DnDs),
                median(Data[Data$Order == 'Lagomorpha',]$DnDs),median(Data[Data$Order == 'Chiroptera',]$DnDs),median(Data[Data$Order == 'Carnivora',]$DnDs),
                median(Data[Data$Order == 'Cetartiodactyla',]$DnDs),median(Data[Data$Order == 'Cetartiodactyla',]$DnDs))
median_KnKs = c(median(Data[Data$Order == 'Eulipotyphla',]$KnKs),median(Data[Data$Order == 'Rodentia',]$KnKs),median(Data[Data$Order == 'Didelphimorphia',]$KnKs),
                median(Data[Data$Order == 'Lagomorpha',]$KnKs),median(Data[Data$Order == 'Chiroptera',]$KnKs),median(Data[Data$Order == 'Carnivora',]$KnKs),
                median(Data[Data$Order == 'Cetartiodactyla',]$KnKs),median(Data[Data$Order == 'Cetartiodactyla',]$KnKs))

median_GenerationLength_d = c(median(Data[Data$Order == 'Eulipotyphla',]$GenerationLength_d),median(Data[Data$Order == 'Rodentia',]$GenerationLength_d),median(Data[Data$Order == 'Didelphimorphia',]$GenerationLength_d),
                              median(Data[Data$Order == 'Lagomorpha',]$GenerationLength_d),median(Data[Data$Order == 'Chiroptera',]$GenerationLength_d),median(Data[Data$Order == 'Carnivora',]$GenerationLength_d),
                              median(Data[Data$Order == 'Cetartiodactyla',]$GenerationLength_d),median(Data[Data$Order == 'Cetartiodactyla',]$GenerationLength_d))
median_AverageGrantham = c(median(Data[Data$Order == 'Eulipotyphla',]$AverageGrantham),median(Data[Data$Order == 'Rodentia',]$AverageGrantham),median(Data[Data$Order == 'Didelphimorphia',]$AverageGrantham),
                           median(Data[Data$Order == 'Lagomorpha',]$AverageGrantham),median(Data[Data$Order == 'Chiroptera',]$AverageGrantham),median(Data[Data$Order == 'Carnivora',]$AverageGrantham),
                           median(Data[Data$Order == 'Cetartiodactyla',]$AverageGrantham),median(Data[Data$Order == 'Cetartiodactyla',]$AverageGrantham))
number_of_species = c(length(unique(Data[Data$Order == 'Eulipotyphla',]$Species)),length(unique(Data[Data$Order == 'Rodentia',]$Species)),length(unique(Data[Data$Order == 'Didelphimorphia',]$Species)),
                      length(unique(Data[Data$Order == 'Lagomorpha',]$Species)),length(unique(Data[Data$Order == 'Chiroptera',]$Species)),length(unique(Data[Data$Order == 'Carnivora',]$Species)),
                      length(unique(Data[Data$Order == 'Cetartiodactyla',]$Species)),length(unique(Data[Data$Order == 'Cetartiodactyla',]$Species)))

Genetics_Ecology = data.frame(Order,median_GenerationLength_d,median_AverageGrantham,median_DnDs,median_KnKs,number_of_species) 
Genetics_Ecology = Genetics_Ecology[order(Genetics_Ecology$median_GenerationLength_d),]

All_statistics = c(GenerationLength_FractionOfSyn_fit$estimate,GenerationLength_FractionOfNonsyn_fit$estimate,GenerationLength_SummOfAllGrantham_fit$estimate,
                   GenerationLength_MedianOfAllNonsyn_fit$estimate,GenerationLength_AverageGrantham_fit$estimate,GenerationLength_KnKs_fit$estimate,GenerationLength_DnDs_fit$estimate)

All_pvalue = c(GenerationLength_FractionOfSyn_fit$p.value,GenerationLength_FractionOfNonsyn_fit$p.value,GenerationLength_SummOfAllGrantham_fit$p.value,
               GenerationLength_MedianOfAllNonsyn_fit$p.value,GenerationLength_AverageGrantham_fit$p.value,GenerationLength_KnKs_fit$p.value,GenerationLength_DnDs_fit$p.value)

Variable = c('FractionOfSyn','FractionOfNonsyn','SummOfAllGrantham','MedianOfAllNonsyn','AverageGrantham','KnKs','DnDs')

Corr_with_GenLen = data.frame(Variable, All_statistics,All_pvalue) 

boxplot(Data[Data$Order == "Eulipotyphla",]$KnKs,+
          Data[Data$Order == "Rodentia",]$KnKs, +
          Data[Data$Order == "Chiroptera",]$KnKs, +
          Data[Data$Order == "Carnivora",]$KnKs, +
          Data[Data$Order == "Primates",]$KnKs, +
          Data[Data$Order =="Cetartiodactyla",]$KnKs, names = c('Eulipotyphla',"Rodentia",'Chiroptera','Carnivora',"Primates",'Cetartiodactyla'), outline = FALSE, xlab = "GenLenght", ylab = "KnKs",notch = TRUE)

boxplot(Data[Data$Order == "Eulipotyphla",]$AverageGrantham,+
          Data[Data$Order == "Rodentia",]$AverageGrantham, +
          Data[Data$Order == "Chiroptera",]$AverageGrantham, +
          Data[Data$Order == "Carnivora",]$AverageGrantham, +
          Data[Data$Order == "Primates",]$AverageGrantham, +
          Data[Data$Order =="Cetartiodactyla",]$AverageGrantham, names = c('Eulipotyphla',"Rodentia",'Chiroptera','Carnivora',"Primates",'Cetartiodactyla'), outline = FALSE, xlab = "GenLenght", ylab = "AverageGrantham",notch = TRUE)

TempData = subset(Data, Data$Order == 'Eulipotyphla'|Data$Order =="Rodentia"|Data$Order =='Chiroptera'|Data$Order =='Carnivora'|Data$Order =="Primates"|Data$Order =='Cetartiodactyla')

ggplot(TempData, aes(x=Order, y=KnKs,fill=Order),notch=T) + 
  geom_violin(trim=FALSE)+ 
  stat_summary(fun=median, geom="crossbar", size=0.2, color="red")+
  labs(title="KnKs~GenLen",x="Order", y = "KnKs")+
  theme_classic()

ggplot(TempData, aes(x=Order, y=AverageGrantham,fill=Order),notch=T) + 
  geom_violin(trim=FALSE)+ 
  stat_summary(fun=median, geom="crossbar", size=0.2, color="red")+
  labs(title="AverageGrantham~GenLen",x="Order", y = "AverageGrantham")+
  theme_classic()

median(Data[Data$Order == "Primates",]$GenerationLength_d)

## [1] 3820.017

median(Data[Data$Order =="Cetartiodactyla",]$GenerationLength_d)# продолжительность жизни китопарнокопытных чуть больше, чем у приматов

## [1] 3595.33

# Carnivora Cetartiodactyla Chiroptera Eulipotyphla Primates Rodentia 

nrow(Mammalia)

## [1] 348

for (i in 1:nrow(Mammalia))   {  Mammalia$FamilyShort[i] = paste(unlist(strsplit(Mammalia$family[i],''))[c(1:3)],collapse = '')  }
for (i in 1:nrow(Mammalia))   {  Mammalia$OrderShort[i] = paste(unlist(strsplit(Mammalia$Order[i],''))[c(1:3)],collapse = '')  }

library(boxplotdbl) # install.packages('boxplotdbl') #здесь надо разбиение по КК 0 и 1
X = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$AverageGrantham)
Y = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X,ylim = c(0,100), xlim = c(0,20000), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE,  draw.legend = TRUE

X = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$DnDs)
Y = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,15000), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE,  draw.legend = TRUE

X = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$KnKs)
Y = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X,ylim = c(0,1), xlim = c(0,18500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE,  draw.legend = TRUE

X = data.frame(as.factor(Mammalia$OrderShort),Mammalia$AverageGrantham)
Y = data.frame(as.factor(Mammalia$OrderShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,6500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE,  draw.legend = TRUE

X = data.frame(as.factor(Mammalia$OrderShort),Mammalia$DnDs)
Y = data.frame(as.factor(Mammalia$OrderShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,6500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE,  draw.legend = TRUE

X = data.frame(as.factor(Mammalia$OrderShort),Mammalia$KnKs)
Y = data.frame(as.factor(Mammalia$OrderShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,6500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE,  draw.legend = TRUE

##########################